其他
使用 R 语言绘制中国乡镇地图及根据经纬度判断所处的乡镇
为了让大家更好的理解本文内容，欢迎各位培训班会员参加明晚 8 点 的直播课：「使用 R 语言绘制中国乡镇地图及根据经纬度判断所处的乡镇」
在之前的课程中，我给大家分享过使用 R 语言绘制省市区县地图的方法：
使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）：https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c
另外也给大家分享过使用 R 语言进行地理编码以及根据经纬度判断所处省市区县的方法：
使用 R 语言进行地理编码：地址解析经纬度、坐标转换 & 根据经纬度判断所处的省市区县：https://rstata.duanshu.com/#/brief/course/6f8633b486ed49e398c486d0d9f0f59a
乡镇地图的绘制
今天我们再来补充讲解下如何绘制乡镇地图以及如何根据经纬度判断所处的乡镇。
附件中的 town4326 文件夹是 WGS84 坐标系下的乡镇地图数据；town_mini 文件夹是经过编辑后的小地图版本的乡镇地图数据；town_long 则是长版本的。
town_mini 和 town_long 都是“+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs”坐标系的，也就是和之前的地图数据一样。
首先加载相关 R 包并读取数据：
library(tidyverse)
library(sf)
# 中国地图通常使用这样的坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
# 小地图版本
# 读取编辑后的乡镇地图数据
read_sf("town_mini/town_mini.shp") -> townmini
# 线条数据可以使用之前课程中提供的省市区县的：
read_sf("chinaprov2021mini/chinaprov2021mini_line.shp") %>%
filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
select(class) -> provlinemap
# 读取绘图数据
haven::read_dta("2019年各乡镇平均人口密度.dta") %>%
rename(townname = 乡镇, towncode = 乡镇代码) -> df
df
#> # A tibble: 43,366 × 3
#> townname towncode mean
#> <chr> <chr> <dbl>
#> 1 一元街办事处 420102003000 9076.
#> 2 一六镇 431022108000 295.
#> 3 一六镇 440232103000 89.2
#> 4 一农场 130209401000 346.
#> 5 一品街道 500113008000 557.
#> 6 一市镇 330226104000 426.
#> 7 一平垣乡 141002203000 677.
#> 8 一平浪镇 532331105000 72.9
#> 9 一座营镇 220581109000 252.
#> 10 一心乡 230624200000 45.7
#> # ℹ 43,356 more rows
然后就可以绘图了，这部分细节较多，建议结合视频讲解学习：
library(ggspatial)
library(ggnewscale)
townmini %>%
left_join(df) %>%
mutate(mean = if_else(is.na(mean), -1, mean)) %>%
mutate(group = cut(mean, breaks = c(-1, -0.1, 68, 125, 197, 312,
458, 663, 1360, 50000),
include.lowest = T,
labels = c("无数据", "<=68",
"68~125", "125~197",
"197~312", "312~458",
"458~663", "663~1360",
">1360"))) -> mapdata
mapdata %>%
ggplot() +
geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA",
"小地图框格" = "black",
"省份" = "black")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3,
"小地图框格" = 0.3,
"省份" = 0.3)
) +
scico::scale_fill_scico_d(palette = "bamako",
name = "人口密度(人/m2)",
direction = -1) +
annotation_scale(location = "bl",
width_hint = 0.3,
text_family = cnfont) +
labs(title = "2019 年各乡镇平均人口密度",
subtitle = "数据爬取&绘制：微信公众号 RStata",
caption = "数据来源：中国科学院资源环境科学与数据中心") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
annotation_north_arrow(
location = "tr",
which_north = "false",
pad_y = unit(0.1, "cm"),
style = north_arrow_fancy_orienteering(
text_family = cnfont
)
) -> p1
ggsave("pic1.png", width = 10, height = 8.5, device = png)
长版地图的绘制代码类似，数据换成长版的即可：
# 长版
read_sf("town_long/town_long.shp") -> townlong
townlong %>%
left_join(df) %>%
mutate(mean = if_else(is.na(mean), -1, mean)) %>%
mutate(group = cut(mean, breaks = c(-1, -0.1, 68, 125, 197, 312,
458, 663, 1360, 50000),
include.lowest = T,
labels = c("无数据", "<=68",
"68~125", "125~197",
"197~312", "312~458",
"458~663", "663~1360",
">1360"))) -> mapdata
# 线条
read_sf("chinaprov2021long/chinaprov2021long_line.shp") %>%
filter(class %in% c("九段线", "海岸线")) %>%
select(class) -> provlinemap
ggplot(mapdata) +
geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3)
) +
scico::scale_fill_scico_d(palette = "bamako",
name = "人口密度(人/m2)",
direction = -1) +
annotation_scale(location = "bl",
width_hint = 0.3,
text_family = cnfont) +
coord_sf(xlim = c(-2625585.8, 2206964.7) * 1.2) +
labs(title = "2019 年各乡镇平均人口密度",
subtitle = "数据爬取&绘制：微信公众号 RStata",
caption = "数据来源：中国科学院资源环境科学与数据中心") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
annotation_north_arrow(
location = "tr",
which_north = "false",
pad_y = unit(0.1, "cm"),
style = north_arrow_fancy_orienteering(
text_family = cnfont
)
) -> p2
ggsave("pic2.png", width = 8, height = 7, device = png)
根据经纬度判断所处的乡镇
首先读取乡镇矢量数据：
read_sf("town4326/town4326.shp") -> town
工企样本.dta 是大概 1 万条工企数据的样本，含有经纬度变量，使用下面的代码就可以判断每个坐标点所处的乡镇了：
haven::read_dta("工企样本.dta") %>%
filter(!is.na(经度)) %>%
st_as_sf(coords = c("经度", "纬度"), crs = 4326, remove = F) -> dfres
如果数据量不大，可以直接用下面的代码判断：
# dfres %>%
# st_intersection(town) -> dfres
数据量大的时候会非常慢，建议转换成 mycrs 坐标系再运算：
# 建议转换成 mycrs 坐标系再判断：
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
dfres %>%
st_transform(mycrs) %>%
st_intersection(st_transform(town, mycrs)) -> dfres2
dfres2 %>%
st_drop_geometry()
#> # A tibble: 10,218 × 7
#> group gqid 企业名称 经度 纬度 Name code
#> * <dbl> <dbl> <chr> <dbl> <dbl> <chr> <chr>
#> 1 719089 1998169629 天津酿酒厂 117. 39.2 丁字沽街道 1201060040…
#> 2 721426 2008228990 大连水泥集团有限公司 122. 39.3 七顶山街道 2102130210…
#> 3 721426 2009042065 大连水泥集团有限公司 122. 39.3 七顶山街道 2102130210…
#> 4 721426 2014087208 大连水泥集团有限公司 122. 39.3 七顶山街道 2102130210…
#> 5 719068 1998169308 天津市神光新技术开发公司 117. 39.1 万兴街道 1201040050…
#> 6 718486 1999000502 北京金盾印刷厂 116. 39.9 万寿路街道 1101080010…
#> 7 718486 2004003087 北京金盾印刷厂 116. 39.9 万寿路街道 1101080010…
#> 8 718486 2007333827 北京金盾印刷厂 116. 39.9 万寿路街道 1101080010…
#> 9 718486 2003004268 北京金盾印刷厂 116. 39.9 万寿路街道 1101080010…
#> 10 718486 2013140022 北京金盾印刷厂 116. 39.9 万寿路街道 1101080010…
#> # ℹ 10,208 more rows
直播信息
为了让大家更好的理解本文内容，欢迎各位培训班会员参加明晚 8 点的直播课：「使用 R 语言绘制中国乡镇地图及根据经纬度判断所处的乡镇」
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/eea9db6b6cd3420e85ab4d8dbc9b4eda